set tcl_precision 8

puts " "
puts "======================================================="
puts "=          Interaccion magnetica Confinado            ="
puts "======================================================="
puts " "

puts "[code_info]"
set pi 3.14159265359

#############################################################
# ENTRADAS DE LA SIMULACION				    #
#############################################################
#------------------------------------------------------------
#Ratio lbox_z/lbox_x (Se considera lbox_x=lbox_y)	    -
#------------------------------------------------------------
set ratio 0.5

#------------------------------------------------------------
#Fraccion en volumen					    -
#------------------------------------------------------------
set phi 0.05

#------------------------------------------------------------
#Intensidad del campo magnetico				    -
#------------------------------------------------------------
set u0 [expr 100.0/8.0]
set m [expr sqrt($u0)]
puts "U0: $u0 (m: $m)"

#------------------------------------------------------------
#Coeficiente de friccion de ec. Langevin
#------------------------------------------------------------
set gamma 10.0

#------------------------------------------------------------
#Temperatura
#------------------------------------------------------------
set temp 1.0

#------------------------------------------------------------
#Reiniciar simulacion    				    -
#------------------------------------------------------------
set restart "yes"
set nFileRestart "system_ini.dat"

#------------------------------------------------------------
#Analisis						    -
#------------------------------------------------------------
set setAnalisis "no"

#------------------------------------------------------------
#Sistema confinado					    -
#------------------------------------------------------------
set setConfined "yes"

#------------------------------------------------------------
#Imprimimos las configuraciones				    -
#------------------------------------------------------------
set printConfig "yes"
#------------------------------------------------------------
#Numero de particulas					    -
#------------------------------------------------------------
set npart  0

if { $restart == "no" } {
    set volume 0.0
    set objFile [open "posicion.dat" "r"]
    set datos [read $objFile]
    close $objFile
    set records [split $datos "\n"]
    set linea(0) 0
    set maxdiam 0.0
    
    ## Iterate over the records
    foreach rec $records {
	    ## Split into fields on colons
	    set fields [split $rec " "]
	    set tonto ""
	    set n 0
    
	    foreach pp $fields {
		    if { $pp != ""} {
			    set linea($n) $pp
			    #puts "$n $pp"
			    incr n
		    }
	    }
	    set x $linea(0)
	    set y $linea(1)
	    set z $linea(2)
	    
	    set diam $linea(3)
	    if {$diam > $maxdiam } {
		    set maxdiam $diam
	    }
	    set volume [expr $volume + $diam*$diam*$diam*$pi/6.0]
	    incr npart
    }
    set npart [expr $npart -1]
	    
}

if { $restart == "yes"} {
	set volume 0.0
	set objFile [open "system_ini.dat" "r"]
	set datos [read $objFile]
	close $objFile
	set records [split $datos "\n"]
	set linea(0) 0
	set maxdiam 0.0
	
	## Iterate over the records
	foreach rec $records {
		## Split into fields on colons
		set fields [split $rec " "]
		set tonto ""
		set n 0
	
		foreach pp $fields {
			if { $pp != ""} {
				set linea($n) $pp
				#puts "$n $pp"
				incr n
			}
		}
		set x $linea(0)
		set y $linea(1)
		set z $linea(2)
		
		set vx $linea(3)
		set vy $linea(4)
		set vz $linea(5)
		
		set fx $linea(6)
		set fy $linea(7)
		set fz $linea(8)
		
		set diam $linea(9)
		if {$diam > $maxdiam } {
			set maxdiam $diam
		}
		set volume [expr $volume + $diam*$diam*$diam*$pi/6.0]
		incr npart
	}
	set npart [expr $npart -1]
}

#------------------------------------------------------------
#Volumen de la caja					    -
#------------------------------------------------------------
set volume     [expr $volume/$phi]

#------------------------------------------------------------
#Lado de la caja					    -
#------------------------------------------------------------
set lx  [expr pow($volume/$ratio,1.0/3.0)]
set ly $lx
set lz [expr $ratio*$lx]

puts "Lado de caja, lx: $lx $ly $lz"

#exit

#############################################################
# INICIALIZACION                                            #
#############################################################


if {$setConfined == "yes"} {
	#------------------------------------------------------------
	#Confinamiento en eje z					    -
	#------------------------------------------------------------
	puts "Hay confinamiento (barreras en eje z)"
	setmd periodic 1 1 0
} else {
	#------------------------------------------------------------
	#No Confinamiento 					    -
	#------------------------------------------------------------
	puts "No hay confinamiento"
	setmd periodic 1 1 1
}

#------------------------------------------------------------
#Creamos la caja 					    -
#------------------------------------------------------------
setmd box_l $lx $ly $lz

#------------------------------------------------------------
#Paso de tiempo 					    -
#------------------------------------------------------------
set dt 0.0005
set dth [expr 0.5*$dt]
setmd time_step $dt

#------------------------------------------------------------
#Numero de ciclos 					    -
#------------------------------------------------------------
set n_cycle 100000000

#------------------------------------------------------------
#Numero de pasos por ciclo				    -
#------------------------------------------------------------
set n_steps 1

#------------------------------------------------------------
#Para la lista de vecinos 				    -
#------------------------------------------------------------
setmd skin 0.5

#------------------------------------------------------------
#Simulacion NVT			                            -
#------------------------------------------------------------
integrate set nvt

#------------------------------------------------------------
#Simulacion NVT			                            -
#------------------------------------------------------------
thermostat langevin $temp $gamma

#------------------------------------------------------------
#Archivo de estado de la simulacion                         -
#------------------------------------------------------------
set nFileSystem "system_bk.dat"

#------------------------------------------------------------
#Impresion de resumen de la simulacion                      -
#------------------------------------------------------------
puts [setmd box_l]
puts [setmd time_step]
puts [integrate]
puts [thermostat]

#############################################################
# PREPARAMOS LOS CORRELADORES                               #
#############################################################

if {$setAnalisis =="yes"} {
    set base 1.25
    set npc 32
    set nc [expr 1000000/$npc]
    set k 0
    set kp 0
    set tciclo 0
    
    for {set i 0} {$i<500000} {incr i} {
    	set tempi($i) 0
    	set nss($i) 0
    }
    
    for {set i 0} {$i < $nc} {incr i} {
    	for {set j 0} {$j <= $npc+1} {incr j} {
    		set val [expr $tciclo + pow($base, $j)]
    		set tempi($k) val
    		set nss($kp) [expr int($val)]
    		if {$kp >=1 && $nss($kp) == $nss([expr $kp -1])} {
    			set kp [expr $kp -1]
    		}
    		incr k
    		incr kp
    	}
    	set tciclo $nss([expr $kp -1])
    }
    
    set inext $nss(0)
    set icorr 0
	
}

#############################################################
# INTERACCIONES                                             #
# Aqui introducimos la de volumen excluido		    #
# 1- A la Klingemberg                                       #
# 2- Interaccion esfera blanda (Weeks-Chandler-Andersen)    #
#############################################################

#------------------------------------------------------------
#2- Interaccion Klingemberg (BMHTF potential modificado)
#------------------------------------------------------------
set cut 4.0
#inter 0 0 lennard-jones 1.0 1.0 $cut 0. 0.0
inter 0 0 klingemberg $u0 100.0 $cut
inter 0 1 klingemberg $u0 100.0 [expr $cut]


#############################################################
# CREAMOS LAS PARTICULAS                                    #
#############################################################

if {$restart == "no"} {
    set objFile [open "posicion.dat" "r"]
    set datos [read $objFile]
    close $objFile
    set records [split $datos "\n"]
    set linea(0) 0
    set nnn 0
    ## Iterate over the records
    foreach rec $records {
	## Split into fields on colons    
	set fields [split $rec " "]
	set tonto ""
	set n 0
	foreach pp $fields {
	    if { $pp != ""} {
		    set linea($n) $pp
		    incr n
	    }
	
	    set x $linea(0)
	    set y $linea(1)
	    set z $linea(2)
	    set diam $linea(3)
	    set masa [expr $diam*$diam*$diam]
	    if {$nnn >= $npart} {
		    set n $npart
	    }
	    part $nnn pos $x $y $z dipm 1.0 dip 0 0 [expr $m*$masa] diameter $diam mass $masa type 0
	}   
	incr nnn
    }
    set min [analyze mindist]
    puts "Distancia minima: $min"
}
	
if { $restart == "yes"} {
    set volume 0.0
    set objFile [open "system_ini.dat" "r"]
    set datos [read $objFile]
    close $objFile
    set records [split $datos "\n"]
    set linea(0) 0
    set nnn 0
    
    ## Iterate over the records
    foreach rec $records {
	## Split into fields on colons
	set fields [split $rec " "]
	set tonto ""
	set n 0
	
	foreach pp $fields {
	    if { $pp != ""} {
		set linea($n) $pp
		#puts "$n $pp"
		    incr n
	    }
	}
	set x $linea(0)
	set y $linea(1)
	set z $linea(2)
	
	set vx $linea(3)
	set vy $linea(4)
	set vz $linea(5)
	
	set fx $linea(6)
	set fy $linea(7)
	set fz $linea(8)
	
	set diam $linea(9)
	set masa [expr $diam*$diam*$diam]
	if {$nnn < $npart} { 
	    part $nnn pos $x $y $z v $vx $vy $vz f $fx $fy $fz dipm 1.0 dip 0 0 [expr $m*$masa] diameter $diam mass $masa type 0
	}
	incr nnn
    }
    set min [analyze mindist]
    puts "Distancia minima: $min $npart"
}

#############################################################
# CONFINAMIENTO                                             #
#############################################################
if {$setConfined == "yes"} {
    #########################################################
    #Creamos dos paredes				    #
    #########################################################
    puts "maxxx $maxdiam"
    constraint wall normal 0 0 -1 dist   [expr -1.0*$lz] type 1
    constraint wall normal 0 0 1 dist  [expr 0.0] type 1
}

#############################################################
# INTERACCION MAGNETICA                                     #
# <r_cut> <mesh> <cao> <alpha>                              #
#############################################################

#------------------------------------------------------------
#Presencia de un campo externo                              -
#------------------------------------------------------------
#constraint ext_magn_field 0 0 10.0

#------------------------------------------------------------
#Interaccion magnetica P3M                                  -
#------------------------------------------------------------
#inter magnetic 1.0 p3m 5.0 32 3 0.5
#inter magnetic 1.0 p3m tunev2 accuracy 1.e-2

#------------------------------------------------------------
#Interaccion magnetica vecinos cercanos                     -
#------------------------------------------------------------
inter magnetic 1.0 dawaanr


#############################################################
# CREAMOS ARCHIVOS DE OBSERVABLES	                    #
#############################################################

set obsEner [open "energia.dat" "w"]
set obsCorrelador [open "correlador.dat" "w"]

#############################################################
# BUCLE PRINCIPAL			                    #
#############################################################

set i 0
set ncount 0
set idconfig 0
set k 0
while { $k<$n_cycle } {
    set ncount 0	
    integrate $n_steps
#    if {$k%1000 == 0} {
#	    analyze append
#    }
#    if {$k>0} {
#	    if {$k%1000 ==0} {
#		    set pp [analyze vanhove 0 0 10.0 1]
#		    analyze remove
#	    }
#    }
    #Tiempo al que corresponde la configuracion
    set t [expr $k*$n_steps*$dt]
    
    #####################################################
    # ANALISIS 			                    #
    #####################################################
    #puts "$k $inext"
    if {$setAnalisis =="yes"} {
        if { $k >= $inext } {
        	set name "config.dat"
        	set out [open $name "a"]
        	puts $out "treal $t  i $k npart $npart lbox $lx $ly $lz"
        	#puts $out "StressTensor: $sigma(0,0) $sigma(0,1) $sigma(0,2) $sigma(1,1) $sigma(1,2) $sigma(2,2)"
        	for {set i 0} { $i < $npart } {incr i} {
        		puts $out [part $i print pos v f diameter]
    		}
        	flush $out
        	close $out
        	    
        	#puts "$icorr  $k $name"
        	incr icorr
        	set inext $nss($icorr)
        }
    }
    
    #####################################################
    # IMPRIMIMOS COSAS		              	    #
    #####################################################
    
    if {$k%1000 == 0} {
	    set emag [analyze energy magnetic]
	    set ekin [analyze energy kinetic]
	    set etot [analyze energy total]
	    set erep [analyze energy nonbonded 0 0]
	    #set erep [expr $etot-$emag-$ekin]
	set temp [expr $ekin/((3.0/2.0)*$npart)]
	set hora [clock format [clock seconds] -format "%H:%M:%S"]
	puts "($hora) Time: $t, Etotal: [expr $etot/$npart], Temp: $temp"
	puts $obsEner "$t [expr $emag/$npart]  [expr $ekin/$npart] [expr $etot/$npart]  [expr $erep/$npart]"
	flush $obsEner
	#puts "Emagn: $emag"
	set systembk [open $nFileSystem "w"]
	puts $systembk "treal $t u0 $u0  npart $npart lbox $lx $ly $lz"
	puts $systembk "[expr $emag/$npart] [expr $erep/$npart] [expr $ekin/$npart] [expr $etot/$npart] "
	    for {set i 0} { $i < $npart } {incr i} {
			      puts $systembk [part $i print pos v f diameter]
	    }
	    close $systembk
	    
    #		if { $printConfig == "yes" } {
    #        		set out [open "config.dat" "a"]
    #        		puts $out "treal $t tcorr 0  icorr 0  npart $npart lbox $lx $ly $lz"
    #        		puts $out "StressTensor: 0 0 0 0 0 0"
    #        		for {set i 0} { $i < $npart } {incr i} {
    #        			puts $out [part $i print pos v f diameter]
    #        		}
    #        		close $out
    #		}
    }
    
	incr k
}



#############################################################
# CERRAMOS ARCHIVOS DE OBSERVABLES	                    #
#############################################################
close $obsEner
close $obsCorrelador
